Diagrammatic Monte Carlo 



oo 
O 
O 
(N 

X> 
<D 
P-h 

O 



O 

£ 

•4— > 

a 

C 

o 
o 



> 
m 

(N 

(N 

(N 
O 
OO 

O 

> 

X 



Kris Van Houcke 1 ' 2 , Evgeny Kozik 1,3 , N. Prokof'ev 1 ' 3,4 , and B. Svistunov 1 ^ 4 

1 Department of Physics, University of Massachusetts, Amherst, MA 01003, USA 

2 Universiteit Gent, Vakgroep Subatomaire en Stralingsfysica, 
Proeftuinstraat 86, B-9000 Gent, Belgium 

3 Institut fur Theoretische Physik, ETH Zurich, CH-8093 Zurich, Switzerland 

4 Russian Research Center "Kurchatov Institute", 123182 Moscow, Russia 

Abstract. Diagrammatic Monte Carlo (DiagMC) is a numeric technique that al- 
lows one to calculate quantities specified in terms of diagrammatic expansions, 
the latter being a standard tool of many-body quantum statistics. The sign prob- 
lem that is typically fatal to Monte Carlo approaches, appears to be manageable 
with DiagMC. Starting with a general introduction to the principles of DiagMC, 
we present a detailed description of the DiagMC scheme for interacting fermions 
(Hubbard model), as well as the first illustrative results for the equations of state. 

1 Introduction. General Principles 

Diagrammatic expansion (Feynman diagrams) is a powerful generic tool of 
quantum statistics [1]. Mathematically, diagrammatic expansion — for some 
relevant quantity, Q, usually a Green's function — is a series of integrals with 
an ever increasing number of integration variables, 



Here y is a set of parameters which the quantity Q can depend on, £ m indexes 
different terms of the same order m (the term m = is understood as a func- 
tion of y only), and the x's are the integration variables. Structures similar 
to Eq. (1) originate from perturbative expansions of the quantum-statistical 
averages, in which case the functions T> can be represented by Feynman di- 
agrams/graphs, with the graph lines standing for either Green's functions 
(propagators) G or interaction potentials U, so that the whole diagram en- 
codes a certain product of G"s and IPs. It will be essential for our purposes 
to approach the integrations in Eq. (1) in the same way as the diagram order 
m and its topology when defining the notion of a diagram, in contrast to 
analytic treatments where integrations are included into the definition of T>. 

Diagrammatic Monte Carlo (DiagMC) [2] is a technique that allows one 
to simulate quantities specified in terms of a diagrammatic series. In a broad 
sense, it is a set of simple generic prescriptions for organizing a systematic- 
error-free Metropolis-Rosenbluth- Teller type process that samples the series 
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(1) without explicitly performing integrations over the internal variables in 
each particular term. The rules are as follows. 

The function Q{y) is interpreted as a distribution function for the vari- 
able (s) y. The statistical interpretation of Eq. (1) suggests calculating Q(y) 
by a Markov-chain process which samples diagrams stochastically. The value 
of V plays the role of the statistical weight of the diagram (i.e., the prob- 
ability with which the diagram is generated). More precisely, to deal with 
sign- alternating series, one writes T> as a product of its absolute value |2?|, 
which determines the statistical weight of the diagram, and the sign func- 
tion sx> = ±1, which contributes to the statistics for Q each time the term 
T> is generated. Since contributions of different m-th order diagrams are of 
the same order of magnitude and their number grows factorially with m, the 
combined m-th order contribution to Q is the result of nearly complete can- 
cellation of the diagrams of different sign. This is the notorious sign problem 
(see, e.g., [3]). In our case, it causes an exponential scaling of the computation 
time with the maximum diagram order, making it impossible to obtain sen- 
sible results at large m. At this point one faces the problem of extrapolating 
the answer to the m — > oo limit. To make matters worse, the series may be 
asymptotic/divergent. For this reason, the use of resummation techniques [4] 
is a crucial ingredient of DiagMC (see Sees. 4, 5). 

Though DiagMC has to confront the sign problem, just like other exact 
techniques (e.g., the dynamical cluster algorithms [5]), it has an important 
advantage: it works immediately in the thermodynamic limit, and is not sub- 
ject to the exponential scaling of the computational complexity with the sys- 
tem or cluster volume. DiagMC also has a potential for vast improvements in 
efficiency by using standard tricks of the analytic diagrammatic approach, de- 
veloped to reduce the number of diagrams calculated explicitly term-by-term 
by self-consistcntly taking into account chains of repeating parts as, e.g., in 
the Dyson equation [1] . These ideas are the essence of the Bold Monte Carlo 
scheme introduced in Ref. [6] and successfully applied to the Fermi-polaron 
problem in Rcf. [4]. In Sec. 4, we describe the basic steps in this direction for 
DiagMC, such as the use of bold propagators. 

The stochastic sampling of Q(y) consists of a number of elementary sub- 
processes (or updates) falling into two qualitatively different classes: (I) up- 
dates which do not change the number of continuous variables (they change 
the values of variables in T>, but not the form of the function itself), and (II) 
updates which change the structure of V. The processes of class I are rather 
straightforward, being identical to those of simulating a given continuous 
distribution |D(£ m , y, x\, . . . ,x m )\ of the variables x\, . . . ,x m . 

The crucial part is played by type-II updates, arranged to form comple- 
mentary pairs with the detailed balance condition [7] satisfied in each pair. Let 
A transform the diagram £>(£m, V, x\,..., x m ) into V(£ m+n , y, x\, . . . , x m+n ), 
and, correspondingly, its counterpart B performs the inverse transformation. 
For the new variables we introduce a vector notation: x = {x m+ i, . . . , x m+n }. 
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The first step in A consists of selecting some type of diagram transformation 
and proposing x which is generated from a certain probability distribution 
W(x). The form of W(x) is arbitrary but to render the algorithm most ef- 
ficient, it is desirable that W(x) is (i) of a simple analytic form with known 
normalization, and (ii) close to the actual statistics of x in the final diagram. 
The second step consists of accepting the proposal with probability, P a dd{x). 
In B one simply accepts the proposal for removing the variables x with the 
probability P rom (a;). The pair of complementary sub-processes is balanced if 
the following equality is fulfilled [2]: 

P add (:r) = mm{R(x)/W(x), 1} , (2) 

P Icm (x) = mm{W(x)/R(x), 1} , (3) 

where 

R(x) = {pB/p A )\T>{^ m+n ,y,x 1 ,...,x m ,x)/T>(^ m ,y,xi,...,x m )\ , (4) 

with p^ and ps being the probabilities of addressing the sub-processes A and 
B. [Often, pj[ ^ pg is a natural choice.] The above protocol is a straightfor- 
ward generalization of the standard approach [7,8] to the sampling of func- 
tions with a variable number of continuous arguments. 

Generally speaking, the set of DiagMC updates is specific for a given type 
of the diagrammatic expansion, being sensitive to both the topology and the 
representation (say, momentum or coordinate) of the diagrams. The updating 
strategy can be rendered more/less sophisticated, depending on what is being 
optimized: the simplicity or the efficiency. Examples of particular DiagMC 
schemes can be found in Refs. [2,4,9-11]. 

In this work, we consider an interacting many-body Fermi system and 
work with diagrams in the imaginary-time - momentum representation. We 
employ a sophisticated (at first glance) updating scheme based on the worm 
algorithm idea [12] when updating flexibility and efficiency are achieved by ex- 
tending the configurational space (in our case, the space of allowed diagrams). 
As far as we know, the algorithm being presented is the first application of 
the DiagMC method to connected many-body Feynman diagrams. 

Perturbative expansions in the interaction are often divergent. Dyson gave 
a simple physical argument leading to a sufficient condition for an expansion 
scries to diverge: if changing the sign (or phase) of some parameter 7 im- 
plies an abrupt change of the physical state (e.g., the change of symmetry, 
or an instability), then 7 = is the point of non-analyticity and the expan- 
sion in powers of 7 is a priori divergent [13]. In particular, this means the 
diagrammatic expansion with respect to the interaction in continuous space 
is divergent for both bosons and fermions, since both systems collapse at 
negative (on average) interaction. Obviously, the opposite is not necessarily 
true — that is the absence of singularity in the aforementioned sense does not 
guarantee analyticity (and thus convergence at small enough 7) — but a dif- 
ferent reason for a series to be non-analytic (and thus divergent at arbitrarily 



4 



Kris Van Houcke et al. 



small 7) should be rather exotic. Therefore, we shall rely on the physical 
picture to predict the analytic properties of the series at 7 — > 0. 

However, the real power of the analytic diagrammatic technique lies in 
the possibility of reducing infinite sums of diagrams (no matter convergent or 
divergent) to simple integral equations, e.g., the Dyson equation. In addition 
to that, for divergent asymptotic series one can employ generalized resumma- 
tion schemes to obtain the answer outside of the series convergence radius. 
Finally, when the resummation trick fails, the series can be analyzed using 
approximate (biased) methods, which assume a certain function behind the 
series. An example of the successful application of the resummation method 
to a divergent sign alternating scries for polarons can be found in Ref. [4] . In 
this work we demonstrate that the DiagMC method is a viable approach to 
study interacting many-body systems. Even moderate success in this direc- 
tion is extremely important since other Monte Carlo approaches face a severe 
sign problem before the fully controlled extrapolation to the thermodynamic 
limit can be done. 

2 Model and Diagrams 

In what follows we discuss the Fermi-Hubbard model, 

H = J2(s k - fi^a^a^ + U^2n u n lr . (5) 

k,<r r 

Here a^ k is the fermion creation operator with the quasi-momentum k lying 
in the first Brillouin zone, n ar — a) ar a ar , a —\ , \ is the spin projection, and r 
is the integer radius vector on the d-dimensional lattice. The lattice dispersion 
is given by £k = — 2i Yla=i cos(k Q a), with a and t being the lattice spacing 
and the hopping amplitude, respectively; /i CT is the chemical potential for the 
component <r, and U is the on-site interaction strength. Units are chosen such 
that a and t are equal to unity. 



Fig. 1. The diagram structural element. The dashed 
vertical line represents the interaction potential. The 
upper (lower) solid lines represent spin-up (spin-down) 
fcrmionic propagators. 

We utilize the standard Matsubara diagrammatic technique [1]. The dia- 
grams consist of (vertical) dashed lines standing for the pair interaction po- 
tential U and solid lines representing particle propagators (see Fig. 1). 
We adopt a convention that the imaginary-time axis is horizontal, and is di- 
rected from right to left. Since in our case there are no interactions within 
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one and the same component, we require that the upper (lower) end of each 
dashed line corresponds to the spin- up (spin-down) component. With this 
rule, the spin-up lines are distinguished from spin-down ones without explicit 
labeling. 

A free fermionic propagator G^(k, t = r 2 — ti) associated with each 
particle line running from t\ to r 2 is defined by 

e<^)={-c::t ( W"' k,)> t;:- <•> 

I +e ^ k n)£ , if r < , 

with = [l + e /3 ( £k_M<T ^] 1 being the occupation of the state (cr, k) at 
inverse temperature (3 for free fermions on the lattice. A Green's function 
with equal time variables (a closed fermion loop) is understood as G^ (k, r = 
— 0). To obtain the right weight and sign, one also has to ascribe the factor 
(—1) N+Nl / (2w) Nd to each diagram of order N, with Ni being the number of 
closed fermion loops and d the dimensionality of the problem. 

All the physical information we will need is contained in the self-energy 
S a [1] . In analytic treatments, it is convenient to have S a in the momentum- 
imaginary-frequency representation, so that the Green's function is obtained 
from the Dyson equation by simple algebra: 

[gAp^T 1 = [g^Hp,^ 1 -z*{p,0 • (7) 

Numerically, we find it more appropriate to work in the momentum-imaginary- 
time representation, to avoid dealing with poles. This does not create a prob- 
lem with finding G cr (p,^), since E a (p,£) is readily obtained from S a (p,r) 
by a (fast) Fourier transform. 

The class of diagrams contributing to ^(p, r) is defined by the following 
requirements: (i) Each diagram has two special vertices having one open 
fermionic end. These vertices are separated by the time interval r and their 
open ends have the same spin projection a and momentum p which enters 
the diagram at one open end and exits at the other, (ii) Each diagram is 
connected, (iii) Each diagram is irreducible, i.e., it can not be split into two 
disconnected parts by cutting only a single fermionic line. 



3 Updates 

For the sake of algorithmic elegance, we choose to work with closed-loop di- 
agrams that (formally) have no free ends. To collect statistics for E in this 
representation, we fix one propagator, say going from vertex 1 to vertex 2, 
and set its value to 1. Below, we refer to this special propagator as 'measur- 
ing' propagator. Then, if the propagator's momentum is p and the vertices 
have times ti,t 2 respectively, we are measuring E(p,n — r 2 ). Note that 
S(p, t) is an anti-periodic function of r (with r = n — r 2 ), so we can restrict 
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ourselves to collecting statistics for < r < (3. For a given diagram, chang- 
ing the measuring propagator — without changing the diagram structure — is 
straightforwardly done by the Swap update. 

Swap. The update swaps the measuring propagator with one of the regular 
propagators. The new measuring propagator is chosen at random, and its 
value is set to 1. Correspondingly, the physical value of the old measuring 
propagator is restored, so that it becomes just a regular propagator. 

Assume that the initial measuring propagator has to be replaced with 
Ga\p,r), after we swap it with the propagator G^(p',r'). The acceptance 
ratio is given by 

P swap = |G(°>(p,T)/GgV,T')|. (8) 

The problem of generating diagrams for S consists of two main tasks: one 
should be able to (i) change the structure of the diagrams, and (ii) make sure 
that the momentum conservation is satisfied in every vertex. We develop 
a scheme that fulfills these tasks by performing only local updates of the 
diagram. The idea — in the spirit of the worm algorithm — is to introduce 
non-physical diagrams in which momentum conservation is violated in some 
special vertices (called worms) . The smallest required number of such vertices 
is two. All updates changing the topology of the diagram involve worms. 
When worms are deleted from the diagram we return to the physical subspace. 

To be specific, we introduce the following notation. One of the three- 
point vertices in which momentum conservation is violated, if any, is labeled 
X. There is an excess momentum S associated with X, which is defined as the 
difference between the incoming and the outgoing momenta in this vertex (see 
Fig. 2). If present, X always has one, and only one, conjugate vertex M with 
an excess momentum —d, as shown in Fig. 2. Clearly, the distinction between 
X and M. is merely conventional, since replacing 8 with — S interchanges X 
and M. There is an important property of the two worms: once some path 
from I to Ai is known, one can simply remove the worms by propagating the 
excess momentum S along the path. 




p-(p,+qi)=8 p'-GV-m')=-8 

Fig. 2. Worms. Graphically, we can picture the worms as vertices with conserving 
momentum, by connecting the two with an extra line (thread) that transports the 
momentum S from I to M. The indistinguishability of T and M is seen from the 
fact that swapping the labels T and M, and replacing 6^—8 yields the same 
original diagram. 
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We proceed now with the description of the updates. The Create/Delete 
pair switches between the physical and non-physical (worm) diagrammatic 
spaces by creating/deleting I, M. Its role is also to update the diagram's 
momenta. The Add/Remove pair changes the diagram order. In these up- 
dates, we add (remove) a vertex and delete (create) X, M. on different spin 
components at the same time. The self-complementary updates Move and Re- 
connect are responsible for moving the worms in the diagram and changing 
its topology, respectively. 

Create. This update is possible only — being rejected otherwise — when we 
are in the physical sector (I7-sector), that is when X, M are absent. We in- 
troduce X, M. by selecting a propagator (line) at random and adding a mo- 
mentum 5 to it (see Fig. 3). [Similarly, X and M can be also created on the 
two vertices of an interaction line.] 

Creating two worms takes the diagram from the I7-sector to the non- 
physical IF-sector. We define the JF-diagram weight by exactly the same 
rules as for physical ones up to an arbitrary numeric factor, so that the 
acceptance ratio for Create is 

P croatc - 2NC N G^(p + 8 1 r)/{2n) d W(8) G^(p,r) , (9) 

where 2N is the total number of propagators (including the measuring prop- 
agator) and W(S) is the normalized distribution from which 5 is drawn. The 
variable Cn is an extra weighing factor which is assigned to a diagram of 
order N whenever the worm is present. To make transitions between the S- 
and ^-sectors efficiently, we choose Cn — C/2N, where C is a constant. 

The simplest choice for W(S) is to assume a uniform distribution in the 
first Brillouin zone, W(S) = l/(2ir) d . Since X and A4 are indistinguishable, 
we find it convenient to require that the leftmost worm is labeled X, and the 
rightmost one M. (see Fig. 3). 



P 2 P Pi .. P 2 ' P ™ Pi 

— A — i 4 1 < p=p+8 — «-Q— <^-«— 

q 2 | jqi q 2 | fqi 

i i i i 

P=Pi-q, =P 2 +q 2 p=p,-q l +5=p 2 +q 2 +8 
Fig. 3. Create update. 



The rest of the updates (except for Remove) apply only — being automat- 
ically rejected otherwise — when X, M. are present. 

Delete. Here, we first check if X and M are connected by a single line and 
proceed only if they are. We propose to remove X and M by adding — 6 (+5) 
to the momentum of the connecting line when it is incoming (outgoing) for 
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X. If the line connecting X and M is a propagator, the acceptance ratio for 
Delete is given by the inverse of Eq. (9), 

Pdcictc - (2v) d W{8) G(°> (p, r) / 27V C N G<°> (p ± S, r) . (10) 

Care should be taken in the special case when X and Ai are connected 
by two lines (forming a closed fermion loop). A factor of 2 [1/2] should be 
included in Eq. (10) [Eq. (9)], because the excess momentum 8 can be at- 
tributed to any of the two lines with probability 1/2. 




Fig. 4. Example of Move left update. 

Move. We can move both X and M. along the lines through the whole 
diagram. Moving a worm along a line is only allowed if there is no other worm 
on this line. [Otherwise, the update is equivalent to Delete.} For example, to 
move T left (Fig. 4) we must add <5 to the line on the left from T. This will 
restore the momentum conservation in X and create a momentum discrepancy 
8 in the vertex left from X. This is our new X, temporarily denoted as X, to 
avoid confusion with the old one. 

The acceptance ratio of this move is given by 

Pmovc = G(°)(p 2 +<5,T)/G(°>(p 2 ,T) , (11) 

with G(°)(p2,t) the outgoing propagator of I before the update. Moving X 
to the left or right is chosen with equal probability. 

Graphically, we can picture the worms as vertices with conserving momen- 
tum, provided we draw an artificial thread directly connecting the worms and 
transporting the momentum 8 from X to M (see Fig. 2). In this language, 
removing the worms means merging the thread with some path formed by the 
lines connecting the worms and adding the thread momentum to their mo- 
menta. If only a part of the connecting path is merged with the thread, this 
results in moving the worms. This way of thinking visually tells us whether 
we have to add or subtract 8 on lines while moving the worms. 

Moving I, M "up" or "down" along an interaction line, we switch to a 
different spin component. For example, to move X down (up) (see Fig. 5), 
we have to add <5 {—8) to the dashed line (interaction) momentum. The 
acceptance ratio for moving X down is given by 



Pmovc — (7(q 2 + s)/u(<a) 



(12) 
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P4 P 3 
Fig. 5. Example of Move down update. 



|q 2 +8 

—*-0-~4— 

P 4 j P3 



Obviously, the move right /lcft/up/down updates are the same for M with 
the only difference that instead of adding ±8 to the lines we have to subtract 
it. 

Reconnect. This simple update can be performed — and is automatically 
rejected otherwise — only if X and M. occupy the same spin component. One 
swaps the incoming end of X with that of M. It is important, however, that 
the update changes the excess momentum 8 associated with the worms. [It is 
also worth noting that without invoking the worms such an update would be 
generally impossible by momentum conservation.] If before the update the 
momenta on the incoming ends of X and M. were pi and P2 respectively, 
then after the update the new excess momentum is given by 8 = 8 + P2 — pi- 

The acceptance ratio for Reconnect is 

p GWjp^TM-n) G(°)(p 2 ,r J -r 2 ) 

rCC G(0)( Pl ,Tl-T 1 )G(0)(p 2 ,T A ,-T 2 ) 

with n, tm the imaginary times of the vertices at which X and M. are located, 
respectively. The variable n (T2) gives the imaginary time of the vertex from 
which a propagator runs to the vertex of X (M) with momentum pi (P2), 
before the update. 



(13) 



J Pi P,-« Pi 

— «-Q— «■ — « — — «— ■ — * — r- «- 

\ \ i i 



i 6 



11 i : 

sis. 1 

— 4 i— 4 — < — 1 4 

M P2 P2+8 

Fig. 6. Add update. 



Add. This update, adding a new interaction line, is only possible — being 
automatically rejected otherwise — if X and M have different spin compo- 
nents. For definiteness, we insert the new interaction line between the propa- 
gator coming into X and the one coming into M, as in Fig. 6. At a randomly 
chosen time, these lines are broken and the four-point vertex is inserted in 
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the breaks. In this particular realization, there is no freedom in choosing the 
momentum along the new dashed line: we set it equal to 6, which means that 
the update inevitably deletes X and M . 
The acceptance ratio for Add is 



-Padd = 



U(S) 



(JV + 1) Wi(r) C n 



G ( f\p 1 -S,r T -r) G^(pi,r-Ti) 



(0), 



Gt 0) (pi,Ti - n) 



G { °\p2 + S,t m -t) G7(p 2 , t-t 2 ) 



(0), 



G'j 0) (p2,T A1 - T 2 ) 



(14) 



where AT is the order of the diagram and W\ (r) is the distribution from which 
the time r of the new interaction vertex is drawn. We have assumed that I 
occupies the spin-up component, without loss of generality. 

Remove. This update is a straightforward inverse of the previous one. It 
can be performed (accepted) only in the Z'-sector. It simultaneously removes 
an interaction line and creates a pair of worms. If the removed line had 
momentum q, for I and A4 we shall have 8 = q. 

In view of the indistinguishability of 1 and A4, we require that 2 (M) is 
created on the spin-up (spin-down) line. The acceptance ratio of this move 
is given by the inverse of Eq. (14) with N being the final diagram order 
now. Adding or removing an interaction line can also involve a change of the 
measuring propagator when this propagator is the incoming propagator of X 
or M. 



4 Useful Tricks and Relations 

Connectivity and irreducibility. The updates discussed above are con- 
structed in such a way that disconnected physical diagrams are never sam- 
pled. In the JC-sector, the only update that can change the topology of the 
diagram is the Swap move. However, swapping the measuring propagator 
cannot create a disconnected piece. When we are dealing with a diagram 
in the M^-sector, Add-Remove and Reconnect are the only updates that can 
change the topology of the diagram. These updates can create two discon- 
nected pieces in the Vy-sector. In this case, however, each disconnected piece 
will contain one worm end. Going back to the Z'-sector amounts to bringing 
the worm ends together via Add or Reconnect, at the same time restoring the 
connectivity. 

To keep track of the irreducibility of the diagrams, we use a hash table of 
lines momenta. The idea is that whenever the diagram is reducible, there is at 
least one propagator which carries the total momentum of the diagram. The 
rcducibility of the diagram can then be established easily by checking the hash 
table and finding a momentum equal to that of the measuring propagator. 
We also check irreducibility right before each measurement of the self-energy 
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Bold propagators. The purpose of this trick is to reduce the space 
of diagrams sampled by Monte Carlo. As already mentioned, this becomes 
essential when one is dealing with a sign-alternating series, in which the sign 
problem scales exponentially with the number of diagrams. Basically, the idea 
follows from Dyson's equation, which allows one to reconstruct the complete 
function G from its elementary building blocks S . Next, in the series for S, 
one has to eliminate all diagrams already accounted for by the replacement 
q(o) _^ q done for all propagator lines. At this point the scheme becomes self- 
consistent. In addition, one can employ geometrical series to define screened 
interaction lines, and, correspondingly, eliminate all diagrams with simple 
fermionic loops. At the time of writing we are in the process of implementing 
the full version of the bold-line trick. 

In the simplest implementation, we build the diagrams on a modified 
Green's function G instead of G^°\ The function G is obtained from G^ 
by incorporating the lowest order "tadpole" diagram (i.e., the mean-field 
solution) self-consistently. This reduces the number of sampled diagrams, 
since all diagrams that can be split into two disconnected parts by cutting 
a single interaction line should be left out, being already taken into account 
in the new Green's function. Mathematically, G is obtained from Eq. (7) by 
simply shifting the chemical potential, [i a — ► /j, a — I]j±\ where — Un a is 
the self-energy in the lowest order and n„ is the density of the component a. 
In principle, updates creating tadpole-contributions can simply be rejected. 
However, to keep the Markov sampling of diagrams ergodic, the presence of 
a limited number of closed propagators G CT (k, t = —0) is allowed, but such 
diagrams are excluded from the statistics of S a . 

Resummation of the diagrams. The use of resummation techniques 
is a crucial ingredient of the diagrammatic Monte Carlo approach [4]. We 
found that the resummation methods effectively reduce the error bars, and 
improve convergence of the Monte Carlo results. In general, for any quantity 
of interest — in our case, self-energy — one constructs partial sums 



defined as sums of all terms up to order A* with the A-th order terms being 
multiplied by the factor F^*\ which has a step-like form as a function of N: 
in the limit of large A* and A <C A* the multiplication factors F approach 
unity while for A — > A* they suppress higher-order contributions in such 
a way that the series J2n=i DnF^*^ becomes convergent. The only other 
requirement is that the crossover region from unity to zero has to increase 
with A*. There are infinitely many ways to construct multiplication factors 
satisfying these conditions. And this yields an important consistency check: 
final results have to be independent of the choice of F. In the Cesaro-Riesz 
summation method we have 




(15) 



N=l 



F^ = [(A* - A + l)/N*] s , (Cesaro-Riesz) . 



(16) 
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Here S > is an arbitrary parameter (5=1 corresponds to the Cesaro 
method). The freedom of choosing the value of Ricsz's exponent S can be 
used to optimize the convergence properties of E(N*). Empirically it was 
found that the factor 

N, 

F (N.) = C (N,) J- exp [-{N. + l) 2 /m(N* - m + 1)] , (17) 

m=N 

where C^ N *^ is such that F^*' = 1, often leads to a faster convergence [4]. 

We proceed as follows. With the series truncated at order AT,, we deter- 
mine the physical quantity of interest (say, number density or energy), and 
then extrapolate its dependence on N* to the — ► oo limit. 

Density and energy. The density is obtained from the Green's function 
through the standard relation 

dk 

' K ■G CT (k,r = -0), (18) 



BZ 



(27T)° 



where the integration is over the first Brillouin zone. Analogously, the kinetic 
(hopping) energy is found as 



£(kin) 



/dk 
jz (2^£kG ff (k,r = -0), (19) 



and the potential (interaction) energy is obtained via 



2i? (p°t) / y = lim T ^ T+0 [ 
Je 



d 

£k + fJ-a 

OT 



G a {k,T-r') ,(20) 



where V is the volume of the system. 



5 Illustrative Results 

To illustrate the application of DiagMC, we simulated the equations of state — 
density and energy as functions of chemical potential and temperature — in 
one (ID) and three (3D) dimensions. Diagrammatically, there is no qualita- 
tive difference between ID and higher-dimensional cases. Meanwhile, in ID, 
where fermions are equivalent to hard-core bosons, we have the advantage 
of comparing DiagMC results with very accurate answers obtained with the 
bosonic worm algorithm [14]. 

By the (discussed in the Introduction) Dyson argument, we expect that at 
any finite temperature the thermodynamic functions arc analytic functions 
of U in a certain vicinity of the point U — 0, because this point is not 
singular in the case of the discrete fermionic system. Correspondingly, the 
expansion in powers of U (that, by dimensional analysis, can be understood 
as the expansion in powers of the dimensionless parameter U/T) is supposed 
to be convergent at small enough U, or, equivalently, large enough T. In 
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ID there are no phase transitions at any finite T, so that one can expect 
that the expansion in powers of U is convergent (or at least resummable) 
down to arbitrarily low temperatures. In higher dimensions, an essential non- 
analyticity of thermodynamic functions appears due to phase transitions. 
Hence, the simplest version of DiagMC described in this work is expected to 
work only at high enough temperatures. 

Note that employing resummation techniques is important not only for 
extending the method to divergent series, or improving the series convergence 
properties, but also for estimating systematic errors due to extrapolation from 
a finite number of terms we are able to calculate. 

The simulation results (plotted versus the inverse of the maximum dia- 
gram order iV*) for ID and 3D and different resummation methods are shown 
in Fig. 7. In ID the horizontal straight lines in both energy and density plots 
mark the exact answers obtained from simulations of two-component bosons 
[14]. The computational effort in 3D was about 4800 CPU-hours. Note that 
the error bars in the final answer are dominated by the systematic — > oo 
extrapolation errors which are estimated from the spread of results for dif- 
ferent extrapolation fits and resummation methods. In particular, for density 
in ID the relative uncertainty due to extrapolation is of the order of 5%, 
suggesting that going to higher is desirable, which is only possible with 
further implementation of the bold-line tricks. 

6 Conclusion 

We have developed a diagrammatic Monte Carlo scheme for a system of 
interacting fermions. For illustrative results, we have obtained equations of 
state (density and energy as functions of chemical potential and temperature) 
for the repulsive Hubbard model. In its current version, the scheme applies 
to the range of temperatures above the critical point. Generalizations of the 
scheme to temperatures below the critical point are possible by using one 
of the following two strategies (as well as their combination), (i) A small 
term explicitly breaking the symmetry (relevant to the critical point) can 
be introduced into the Hamiltonian. The phase transition then becomes a 
crossover, and the convergence/resummability of series is likely to take place 
at any temperature, (ii) Anomalous propagators can be introduced. The latter 
trick naturally involves the bold-line technique which is important on its own, 
even above the critical point, since it reduces the number of diagrams and 
thus alleviates the sign problem. 

Tolerance to the sign problem is the single most important feature of 
DiagMC. One can formulate the sign problem as the impossibility of obtain- 
ing results with small error bars for system sizes which allow a reliable and 
controlled extrapolation to the thermodynamic limit. Within the DiagMC 
approach the thermodynamic limit is obtained for free, while, as demon- 
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Fig. 7. The number density n CT and the total energy density E a /V per spin com- 
ponent in ID [panels (a), (b)] with U = 4.0, /ij^ = —0.5, and T = 0.3, and in 3D 
[panels (c), (d)] with U = 4.0, /xf,| = 1-5, and T = 0.5 plotted as functions of the in- 
verse of the maximum diagram order N t . The exact results in ID, n a = 0.2718 and 
E a /V — —0.3613, were obtained from simulations of two-component bosons [14]. 
The results of different resummation procedures are shown — Cesaro sum (circles), 
Cesaro-Riesz with 5 = 2 (triangles), Cesaro- Riesz with 5 = 4 (inverted triangles), 
and resummation with F^* given by Eq. (17) (diamonds), and bare series with- 
out resummation (squares). The solid lines show the best fits [linear for Cesaro, 
polynomial of order 5 for Cesaro-Riesz, and exponential for Eq. (17)]. 

strated in this work, an extrapolation to the infinite diagram order can be 
done sensibly before the error bars explode. 
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